alpha   = 1-gamma;
rho     = 1-1/psi;
theta   = alpha/rho;

% Solve value function
lnWC        = ones(nS,1);
opt         = optimset('Display','off','TolFun',1e-14);
f_lnWC      = @(lnWC) theta*log(beta) + alpha*mu_C + 1/2*alpha^2*sig_C.^2 + log(P*(exp(lnWC)+1).^theta) - theta*lnWC;
lnWC        = fsolve(f_lnWC,lnWC,opt);
lam         = exp(lnWC);


% Pre-compute some SDF related objects
S           = ((lam'+1)./lam).^(theta-1);
Sbar        = diag(S*P');
s           = log(S./Sbar);
Q           = P.*exp(s);
stat_Q      = Q^1000;
stat_Q      = stat_Q(1,:)';


% Risk-free rate
A           = beta^theta * bsxfun(@times,S,exp(-gamma*mu_C+0.5*gamma^2*sig_C.^2));
bond_f      = NaN(nS,120);
bond_1      = ones(nS,1);
for i=1:120
    bond_1      = (P.*A)*bond_1;
    bond_f(:,i) = bond_1;
end
Rf      = 1./bond_f; % [nS,120]

